%% load sample table
tbl_samples = readtable("../20230926_summary_short_samples.xlsx");

%% load measurement & cell numbers of median cell size
opt = detectImportOptions("tbl_medianArea_Human.csv");
opt.VariableTypes = {'char','double','char','double','double','double','double','double','double','double','double','double','double'};
tbl_h = readtable("tbl_medianArea_Human.csv",opt);
species = repmat("Human",height(tbl_h),1);
tbl_h = addvars(tbl_h,species,'Before','dateCode');

opt = detectImportOptions("tbl_medianArea_RousettusAegyptiacus.csv");
opt.VariableTypes = {'char','double','char','double','double','double','double','double','double','double','double','double','double'};
tbl_ra = readtable("tbl_medianArea_RousettusAegyptiacus.csv",opt);
species = repmat("RousettusAegyptiacus",height(tbl_ra),1);
tbl_ra = addvars(tbl_ra,species,'Before','dateCode');

opt = detectImportOptions("tbl_medianArea_NyctalusNoctula.csv");
opt.VariableTypes = {'char','double','char','double','double','double','double','double','double','double','double','double','double'};
tbl_nn = readtable("tbl_medianArea_NyctalusNoctula.csv",opt);
species = repmat("NyctalusNoctula",height(tbl_nn),1);
tbl_nn = addvars(tbl_nn,species,'Before','dateCode');

tbl_species = [tbl_h; tbl_ra; tbl_nn];
clearvars tbl_h tbl_ra tbl_nn

%% iterate over table: automatic
surfaceArea2volume = nan(height(tbl_species),1);
volume2surfaceArea = nan(height(tbl_species),1);
divCellNumPtr = height(tbl_species)/50;
fprintf([repmat('-',1,50) '\n']);

tic;
for idxTbl = 1:height(tbl_species)
    fprintf(repmat('-',1,round(idxTbl/divCellNumPtr) - round((idxTbl-1)/divCellNumPtr)));
    idxDataset = find((tbl_species.species(idxTbl) == tbl_samples.Species) & (tbl_species.index(idxTbl) == tbl_samples.Index));
    % find file
    if ("" == tbl_samples.Animal(idxDataset))
        str = "";
    else
        str = "*";
    end
    dirFile = dir("../../" + tbl_species.dateCode(idxTbl) + "*" + "/" + str + tbl_samples.Animal(idxDataset) + "*" + tbl_species.tempLevel(idxTbl) + "*" + "/" + "Online" + "/" + "M" + tbl_species.measNum(idxTbl) + "_data_export.mat");
    if length(dirFile) ~= 1
        error("Sample data not found!");
    end
    load(fullfile(dirFile.folder,dirFile.name));
    idxCell = find([sample.measData.cellData.number] == tbl_species.cellNum(idxTbl));
    if sample.measData.cellData(idxCell).data.medianArea - tbl_species.area(idxTbl) > 1e-10
        error("Found cell incorrect!")
    end
    surfaceArea2volume(idxTbl) = median(sample.measData.cellData(idxCell).data.surfaceArea./sample.measData.cellData(idxCell).data.volume);
    volume2surfaceArea(idxTbl) = median(sample.measData.cellData(idxCell).data.volume./sample.measData.cellData(idxCell).data.surfaceArea);
end
fprintf('\n');
toc

tbl_species = addvars(tbl_species,volume2surfaceArea,'After','area');
tbl_species = addvars(tbl_species,surfaceArea2volume,'After','area');

%% categorical
tbl_species.species = categorical(tbl_species.species);
tbl_species.tempLevel = categorical(tbl_species.tempLevel);

%% save table
writetable(tbl_species,"tbl_surfArea2volRatio.csv");

%% statistics
tbl_stats  = cell2table(cell(0,11), 'VariableNames', {'Species', 'TempLevel', 'SpeciesTempLevel', ...
    'Area', 'SurfArea2Vol_mean', 'SurfArea2Vol_std', 'SurfArea2Vol_sem', 'Vol2SurfArea_mean', 'Vol2SurfArea_std', 'Vol2SurfArea_sem', 'NoOfCells'});

levels.species = categories(tbl_species.species);
levels.tempLevel = categories(tbl_species.tempLevel);

for idxSpecies = 1:length(levels.species)
    mask = (tbl_species.species == levels.species{idxSpecies});
    tbl_stats = [tbl_stats; {levels.species{idxSpecies}, "All", levels.species{idxSpecies}, ...
        median(tbl_species.area(mask)),...
        mean(tbl_species.surfaceArea2volume(mask)),std(tbl_species.surfaceArea2volume(mask)),nan,...
        mean(tbl_species.volume2surfaceArea(mask)),std(tbl_species.volume2surfaceArea(mask)),nan,...
        sum(mask)}];
    for idxTempLevel = 1:length(levels.tempLevel)
        mask = (tbl_species.species == levels.species{idxSpecies}) & (tbl_species.tempLevel == levels.tempLevel{idxTempLevel});
        tbl_stats = [tbl_stats; {levels.species{idxSpecies}, levels.tempLevel{idxTempLevel}, levels.species{idxSpecies}+" @ "+levels.tempLevel{idxTempLevel}, ...
            median(tbl_species.area(mask)),...
            mean(tbl_species.surfaceArea2volume(mask)),std(tbl_species.surfaceArea2volume(mask)),nan,...
            mean(tbl_species.volume2surfaceArea(mask)),std(tbl_species.volume2surfaceArea(mask)),nan,...
            sum(mask)}];
    end
end

tbl_stats.SurfArea2Vol_sem = tbl_stats.SurfArea2Vol_std ./ sqrt(tbl_stats.NoOfCells);
tbl_stats.Vol2SurfArea_sem = tbl_stats.Vol2SurfArea_std ./ sqrt(tbl_stats.NoOfCells);
tbl_stats.SpeciesTempLevel = categorical(tbl_stats.SpeciesTempLevel);
writetable(tbl_stats,"tbl_surfArea2volRatio_stats.csv");

%% plot
figure(1);clf;
fPos = get(gcf,'Position'); set(gcf,'Position',[fPos(1:3),500]);
errorbar(tbl_stats.SpeciesTempLevel,tbl_stats.SurfArea2Vol_mean,tbl_stats.SurfArea2Vol_std,'.','LineWidth',1.5,'MarkerSize',15);
xlabel('Species @ temp');
ylabel('Surface area / volume (μm^{-1}), mean ± std');
ylim([.95 2]);
for idxTbl = 1:height(tbl_stats)
    text(idxTbl,.98,"{\itn} = "+tbl_stats.NoOfCells(idxTbl),"Rotation",90);
end

saveas(gcf,"fig_surfArea2volRatio_stats",'fig');
saveas(gcf,"fig_surfArea2volRatio_stats",'png');


figure(2);clf;
fPos = get(gcf,'Position'); set(gcf,'Position',[fPos(1:3),500]);
errorbar(tbl_stats.SpeciesTempLevel,tbl_stats.Vol2SurfArea_mean,tbl_stats.Vol2SurfArea_std,'.','LineWidth',1.5,'MarkerSize',15);
xlabel('Species @ temp');
ylabel('Volume / surface area (μm), mean ± std');
ylim([.4 .9]);
for idxTbl = 1:height(tbl_stats)
    text(idxTbl,.43,"{\itn} = "+tbl_stats.NoOfCells(idxTbl),"Rotation",90);
end

saveas(gcf,"fig_vol2surfAreaRatio_stats",'fig');
saveas(gcf,"fig_vol2surfAreaRatio_stats",'png');


%% former analysis
% %% load data: manual
% 
% % human
% load("../../191108_Bob_HumanBlood1/191108_H01_8nlps_RT/Online/M2_data_export.mat")
% cellIdx = find([sample.measData.cellData.number] == 255);
% % 1.3197
% 
% 
% % rou aeg
% load("../../210812_Faruq_2_FLIbats/210812_Bat1_8nlps_RT/Online/M3_data_export.mat")
% cellIdx = find([sample.measData.cellData.number] == 244);
% % 1.6577
% 
% 
% % nyc noc
% load("../../191015_GeraldDoreenBob_Bats1/191015_Bat1_8nlps_RT/Online/M3_data_export.mat")
% cellIdx = find([sample.measData.cellData.number] == 579);
% % 1.8511
% 
% %% check & output
% % check
% median(sample.measData.cellData(cellIdx).data.area)
% % surface area to volume ratio
% median(sample.measData.cellData(cellIdx).data.surfaceArea./sample.measData.cellData(cellIdx).data.volume)
% 
% 
% %% check parameters
% 
% nexttile
% plot(sample.measData.cellData(cellIdx).data.time,sample.measData.cellData(cellIdx).data.area)
% median(sample.measData.cellData(cellIdx).data.area)
% nexttile
% plot(sample.measData.cellData(cellIdx).data.time,sample.measData.cellData(cellIdx).data.surfaceArea)
% median(sample.measData.cellData(cellIdx).data.surfaceArea)
% nexttile
% plot(sample.measData.cellData(cellIdx).data.time,sample.measData.cellData(cellIdx).data.volume)
% median(sample.measData.cellData(cellIdx).data.volume)
% nexttile
% plot(sample.measData.cellData(cellIdx).data.time,sample.measData.cellData(cellIdx).data.surfaceArea./sample.measData.cellData(cellIdx).data.volume)
% median(sample.measData.cellData(cellIdx).data.surfaceArea./sample.measData.cellData(cellIdx).data.volume)
